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In this work we assess the quahty and performance of several novel dissipative particle dynamics integration 
schemes that have not previously been tested independently. Based on a thorough comparison we identify the 
respective methods of Lowe and Shardlow as particularly promising candidates for future studies of large-scale 
properties of soft matter systems. 
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I. INTRODUCTION 

The mesoscopic phenomena of so-called "soft matter" 
physics |[ll 13 Isj, embracing a diverse range of systems in- 
cluding liquid crystals, colloids, and biomembranes, generally 
involve some form of coupling between different characteris- 
tic time- and length-scales. Computational modeling of such 
multi-scale effects requires new methodology applicable be- 
yond the realm of traditional techniques such as ab initio and 
classical molecular dynamics 14^ |5| (the methods of choice 
in the microscopic regime), and phase field modeling |6] or 
the lattice-Boltzmann method |7] (usually concerned with the 
macroscopic regime). 

Dissipative particle dynamics (DPD) i9l flolfTllll2il is a par- 
ticularly appealing technique in this regard. The "particles" 
of DPD correspond to coarse-grained entities, representing a 
collection of molecules or molecular groups rather than in- 
dividual atoms. Coarse-graining leads to soft pair potentials 
allowing the particles to overlap (Forrest and Suter 1 13|). 

Although coarse-graining might also be considered implicit 
in Brownian and Langevin dynamics simulation, DPD offers 
the explicit advantage of a proper description of hydrody- 
namic modes significant in the physical approach towards a 
system's equilibrium. This is achieved in DPD by implement- 
ing a thermostat in terms of pairwise random and dissipative 
forces such that the total momentum of the system is con- 
served. Due to these reasons, DPD has been used in studies 
covering a wide range of aspects in soft matter systems, in- 
cluding the structure of lipid bilayers ifljlfisll . self-assembly 
lll6ll . and the formation of polymer-surfactant complexes 1 17|. 

In practice, the pairwise coupling of particles through ran- 
dom and dissipative forces makes the integration of the equa- 
tions of motion a nontrivial task. The main difficulty arises 
from the dissipative force, which depends explicitly on the 
relative velocities of the particles, while the velocities in turn 
depend on the dissipative forces. An accurate description of 
the dynamics requires a self-consistent solution. 

The considerable computational load associated with this 
task has motivated the development of schemes |12, 18, 19, 
I20ll2lll23.l23ll2^ which account for the velocity dependence 
of dissipative forces in some approximate manner, allowing 



the integration to be carried out to a sufficient degree of com- 
putational efficiency. The search for a satisfactory such in- 
tegration scheme is ongoing, since many of the recent pro- 
posal have been found to exhibit non-physical behavior, such 
as systematic drift in temperatiire, and artificial structures in 
the radial distribution function j2(Al2lll221 . 

In order to overcome these problems, a number of new in- 
tegration schemes for DPD simulations have been developed 
in the past few years. Self-consistent determination schemes 
exist on the one hand |l20l|2ll|22|], but these are rather elabo- 
rate. 

Alternative proposals include (i) a parameterization of the 
integrator based on the_specific application being modeled by 
den Otter and Clarke |23], (ii) operator splitting by Shardlow 
ll24ll . and (iii) an elegant Monte Carlo-based approach due to 
Lowe 1 25 1 which completely avoids the problems arising from 
random and dissipative forces as it does not use random or 
dissipative forces at all. 

In this article, we apply these schemes respectively to spe- 
cific model systems, with the objective of assessing their rel- 
ative performance. To the authors' knowledge, the latter three 
have yet to be tested and compared independently. The self- 
consistent approach has been tested, although not extensively, 
and the previously tested so-called DPD-VV (DPD version of 
the "velocity- Verlet" scheme) will be used here as a bench- 
mark. By monitoring a number of physical observables in- 
cluding temperature, radial distribution function, radius of 
gyration for polymers, and tracer diffusion, we find that the 
methods by Lowe 125.1 and Shardlow |24] give the best over- 
all performance and are superior also to the integrators tested 
in previous studies As will be discussed in the last 

section, a direct comparison between these two is not straight- 
forward, since they are based on essentially different concep- 
tual views of dissipative particle dynamics. 

The paper is organized as follows. First, we briefly summa- 
rize the dissipative particle dynamics method and introduce 
the three model systems used here. In Sect. |in]we describe 
the integration algorithms and the update schemes in detail, 
and in Sect. II VI we present the results from the tests. Finally, 
in Sect. |V]the findings and their relevance are discussed. 
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II. METHODS AND MODELS 

Below we give a short summary of the dissipative particle 
dynamics method and describe the three model systems used 
in this work. For more thorough accounts on DPD, see e.g. 
Refs. Elllll. 



A. Dissipative Particle Dynamics 

Dissipative particle dynamics describes a system in terms 
of N particles having masses m^, positions fi, and velocities 
Vi. Interactions are composed of pairwise conservative, dissi- 
pative, and random forces exerted on particle i by particle j, 
respectively, and are given by 

pc ^ i^(=Vr - W ■ 

Fjj = -1 (rtj) (vtj ■ e^j) , (1) 

where f^j = fi - fj, = \fi-j\, eij = nj/vij, and 
Vij = Vi — Vj . The variables 7 and a are the strengths of the 
dissipative and random forces, respectively. The are sym- 
metric Gaussian random variables with zero mean and unit 
variance, and are independent for different pairs of particles 
and different times. The condition = is employed to 
ensure momentum conservation. 

(c) 

The pairwise conservative force is not specified by the 
DPD formulation and it can be chosen to include any forces 
that are appropriate for a given system, such as van der Waals 
and electrostatic interactions. In addition, it is important to 
notice that it is completely independent of the random and 
dissipative forces. Since one of the main motivations for using 
DPD is to be able to simulate systems at coarse-grained, or 
mesoscopic, scales, the conservative force is often chosen to 
be soft repulsive. Here, we use the "classical" DPD choice, 
i.e., soft repulsive conservative forces in the cases of model A 
and model B, and hard Lennard- Jones interactions combined 
with harmonic spring forces in the model polymer system. 

In contrast to the conservative force, the random and dissi- 
pative forces are not independent, but are coupled through a 
fluctuation-dissipation relation. This coupling is due to the re- 
quirement that in thermodynamic equilibrium the system must 
have canonical distribution. The necessary conditions were 
first derived by Espanol and Warren in 1995 using a Fokker- 
Planck equation LI 0.1 . 

The requirement of canonical distribution sets two condi- 
tions linking the random and dissipative forces in Eq. Q. 
The first one couples the weight functions through uj^ [rij ) = 
[u)^{rij)Y, and the second one the strengths of the random 
and dissipative forces via — 27fcsT*. The latter condition 
fixes the temperature of the system T* (ks being the Boltz- 
mann constant) and relates it to the two DPD parameters 7 
and <T. 

Like classical molecular dynamics (MD) simulations, DPD 
allows studies of dynamical properties since the time evolu- 
tion of particles can be described by the Newton's equations 



of motion 

dfi ~ Vidt, 

dv, = J-{Ffdt + Ffdt + FfVdi). 

Here Ff = '^i^j Fij is the total conservative force acting 

on particle i, and Ff and Ff are defined correspondingly. It 
is important to notice that the velocity increment due to the 
random force in Eq. (|2j has the factor Vdt instead of dt. It 
can be justified by a Wiener process as in stochastic differen- 
tial equations. Here, it suffices to notice that physically the 
Wiener process models intrinsic (thermal) noise in the sys- 
tem and provides the simplest approach to modeling Brownian 
motion using stochastic processes (see Refs. 1 10, 12] for a de- 
tailed discussion). The above continuous-time version of DPD 
satisfies detailed balance and describes the canonical NVT en- 
semble. In practice, however, the time increments in Eq. Q 
are finite and the equations of motion must be solved by some 
integration procedure. We will return to this issue in Sect.lIIII 



B. An alternative approach to DPD 

Dissipative particle dynamics described above can be 
thought of as a momentum conserving thermostat that allows 
one to study a system within the NVT ensemble with full hy- 
drodynamics. The key features are therefore momentum and 
temperature conservation. As discussed above, momentum 
conservation arises naturally from pairwise forces. Tempera- 
ture conservation, in turn, arises from the random and dissipa- 
tive forces that are chosen to satisfy the fluctuation-dissipation 
theorem. 

An alternative approach was formulated by Lowe in 1999 
It does not use dissipative or random forces at all, yet 
provides the same conservation laws and is similar in spirit 
to DPD as it is aimed for studies of coarse-grained models 
in terms of soft interactions. In Lowe's method, one first in- 
tegrates Newton's equations of motion with a time step At, 
and then thermalizes the system using the Andersen thermo- 
stat | 26j for pairs of particles. We will discuss this method in 
detail in Sect. lmEl 

Lowe's approach is appealing for a number of reasons. First 
of all, since there are no dissipative forces we can assume that 
this method does not suffer from the same drawback as DPD: 
While DPD requires a self-consistent solution of the equa- 
tions of motion, Lowe's approach is easier to use and most 
likely performs well even with integration schemes that are 
commonly used in classical MD simulations. Secondly, the 
rate of how often the particle velocities are thermalized may 
be varied over a wide range, which implies that the dynamical 
properties of the system can be tuned in a controlled fashion. 
Lowe has demonstrated this idea by showing how some di- 
mensionless variables (such as the Schmidt number Sc) can 
be tuned to match values found in actual fluids I25i] . 
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C. Model systems 

In this study, we evaluate the performance of a number of 
novel integration schemes that have been recently suggested 
for large-scale DPD simulations (see Sect.lIIH. We test these 
integrators using three different model systems. The first two 
are based on a 3D model fluid system with a fixed number 
of identical particles. The first of them is aimed to clarify the 
performance of integration schemes in weakly interacting sys- 
tems dominated by the random and dissipative forces, while 
the second model is relevant for systems in which the conser- 
vative interactions are of major importance. Finally, to gain 
insight into problems associated with hybrid models in which 
both soft and hard interactions are included, we consider a 
model of an individual polymer chain in a hydrodynamic sol- 
vent. 



1. Model A 

Model A describes the case characterized by the absence 
of conservative forces (Fj^ — 0). This choice corresponds 
to an ideal gas and it is customarily called "ideal DPD fluid". 
The reason for using this model is that it provides us with 
some exact theoretical results that can be compared to results 
from model simulations. Here, the dynamics of the system 
arises only from thermal noise and dissipative coupling be- 
tween pairs of particles. In DPD simulations, the random 
force strength is chosen to be ct = 3 in units of fc^T*, and 
the strength of the dissipative force 7 is then determined by 
the fluctuation-dissipation relation CT^ — 2jkBT*. 

The random and dissipative forces are chosen to be soft- 
repulsive. 



for nj < rc 
for rij > rc 



(3) 



with a cut-off distance Vc fl^ and uj^{rij) = uj{rij). This 
is the most common choice in DPD simulations and it has 
been used in recent investigations of integration schemes i2ll 
l22il . thus allowing for a comparison of the present results with 
those of previous works. Although Eq. (|3j has been used in 
virtually all published studies using DPD, it should be noted 
that the fluctuation-dissipation theorem does not specify the 
functional form of the weight function. The simple form of 
Eq. (|3} simply provides a convenient choice. 

In our simulations, a 3D simulation box of size 10 x 10 x 10 
with periodic boundary conditions is used. The length scale is 
defined by setting rc = 1, and a particle number density p = 
4. The energy scale is defined by setting the desired thermal 
energy to unity via fc^T* — 1. All particles are identical, and 
thus m, — m for all i. 



we choose to be of the form Fl^\rij) = Abj{rij). The am- 
plitude of the force was chosen Xoht A = 25. This functional 
form for the conservative force is by far the most common 
choice in DPD simulations. 



3. Model polymer system 

The last model system considered in this work describes 
an individual polymer chain in an explicit hydrodynamic sol- 
vent. Our interest in a system of this kind originates from 
the fact that various soft matter systems such as liquid crys- 
tals and lipid bilayers are composed of particles which are es- 
sentially chain-like molecules. DPD serves well for studies 
of these systems due to the hydrodynamic nature of the sol- 
vent which plays an important role in various soft matter sys- 
tems. However, while it is often desirable to describe chain- 
like molecules on a molecular level by hard interactions, com- 
plemented with bending and torsional potentials to account for 
the most relevant microscopic degrees of freedom, the solvent 
can often be described on a simpler level in terms of soft pair 
potentials. Thus, one possibility for efficiently modeling poly- 
meric systems is a hybrid approach of chain-like molecules in 
a coarse-grained solvent. 

The gain of using a hybrid approach in complex macro- 
molecular systems is evident. It allows one to reduce the com- 
putational burden of dealing with an explicit solvent, while the 
molecular description of macromolecules is still accounted for 
in detail. However, the practical implications of including 
both hard and soft interactions in a model are not well un- 
derstood. In a previous study for an ensemble of spherical 
particles described by hard conservative and soft dissipative 
forces, we found certain features which differentiated integra- 
tion schemes from each other [22.1 . However, a full study of 
the performance of integration schemes within a true hybrid 
approach of a macromolecular system has been lacking up till 
now. 

This model system aims to quantify the effects of integra- 
tion schemes under conditions that combine both soft and 
hard interactions for a model polymer system. Here, the idea 
is to optimize the efficiency of the model by using a mini- 
mal approach. We thus describe the solvent as an ensemble 
of identical particles which interact via soft pairwise forces 
and satisfy momentum conservation, while the polymer chain 
is described on a more microscopic level in terms of (hard) 
Lennard- Jones interactions and harmonic springs. 

The linear polymer chain is described as a chain of M 
monomers connected by harmonic bonds whose potential fol- 
lows the form 



2. Model B 



arm 2 



i = 2,3,. 



,M, 



(4) 



Model B is a simple interacting DPD fluid. Its main differ- 
ence to model A is the presence of a conservative force, which 



with a spring constant k = 7. Within the chain, the conser- 
vative monomer-monomer interactions are given by the trun- 



cated and shifted Lennard- Jones potential 




< 



(5) 

such that the potential is purely repulsive and decays smoothly 
to zero at Vc- We choose £ ~ 2^^/^ and e = ksT*, and there- 
fore Tc = £2^/^ = 1. The pairwise conservative force acting 
on a monomer due to other monomers in a chain therefore 
follows directly from F'~' = — V(C/harm + C^lj)- The dissipa- 
tive and random forces acting on the monomers are chosen to 
follow Eqs. Q and with a — 3. 

The monomer- solvent and solvent-solvent interactions are 
described as in model A with cr = 3 (and A ~ 0), i.e., the ran- 
dom and dissipative parts are used as a momentum conserving 
thermostat. 

The justification for this choice of interactions lies in a wish 
to clarify the size of artifacts due to integration schemes in a 
case, where the polymer chain is described on a molecular 
level in terms of hard interactions, while the solvent is coarse 
grained as much as possible and is thus described by an ideal 
gas. The coupling between the solvent and the polymer chain 
comes from the dissipative forces which give rise to hydro- 
dynamic modes. This eventually results in a minimal model 
of a polymer chain with full hydrodynamics under good sol- 
vent conditions. This was confirmed by studying the scaling 
behavior of the radius of gyration. 

We consider polymers of size AI ~ 20 and use a 3D box of 
size 10 X 10 X 10 with periodic boundary conditions, where 
the length scale is defined by Tc = 1. The particle number 
density is chosen to be p = 4 while the energy scale is defined 
by setting the desired thermal energy to unity via ksT* = 1. 
All particles are identical, and thus rrii — m for both solvent 
and monomer particles. 



4. Choice of random numbers 

In the present work, uniformly distributed random numbers 
u e t7(0, 1) are used such that = V3(2m — 1). This ap- 
proach is highly efficient and yields results that are indistin- 
guishable from those generated by Gaussian random numbers 
llT2il . However, in the case of Lowe's approach, the ^[^^ used 
are true Gaussian random numbers. 



III. INTEGRATORS 

The integration schemes tested in this work have been cho- 
sen from the most recent ones that have been suggested in 
the literature but not tested and compared to other methods. 
They complement each other in the sense that the velocity 
dependence of the dissipative forces is accounted for in all 
cases, but the approaches differ substantially. We feel that all 
of the integrators considered here are promising candidates 
for large-scale simulations of soft matter systems. However, 
due to the lack of comparative studies in which all of these 
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PfAt + FfAt + Ff 



(1) V, 

(2) r, 

(3) Calculate {f, }, i^^{rj , Vj }, i^«{f, } 



AtV 



(4a) v° 
(4b) V, 



Vi + ■ 



Ff^At + F/^ 



At] 



^ 2 m * 



(5) Calculate F/' {f J , iT, } 

(6) Calculate physical quantities 



TABLE I: Update scheme for DPD-VV and its self-consistent ver- 
sion SC~VV. In the case of DPD~VV, steps (4b) and (5) ai'e done 
only once during a single time step. For the self-consistent integrator 
SC~VV, the loop over steps (4b) and (5) is repeated until the instan- 
taneous temperature has converged to its limiting value. 



schemes would have been tested on equal footing, their rel- 
ative performance has remained an open question. Here, we 
clarify this situation. 



A. Velocity- Verlet based integration scheme DPD-VV 

In Table U we summarize the simplest integrator tested in 
this study. DPD-VV i2lll23l is based on the standard molec- 
ular dynamics velocity- Verlet algorithm l"?, |^ Esil which 
is a time-reversible and symplectic second-order integration 
scheme. These properties can be proven by a straightfor- 
ward application of the Trotter expansion f4]. The standard 
velocity- Verlet has been shown to be relatively accurate in 
typical MD simulations especially at large time steps |28|]. 
The simplicity and good overall performance of the velocity- 
Verlet algorithm thus makes it a good starting point for further 
development. 

We use the acronym DPD-VV for the modified velocity- 
Verlet. DPD-VV differs from the standard velocity-Verlet 
scheme in one important respect. As discussed above, the 
dissipative forces in DPD depend on the velocities, which 
in turn are governed by the dissipative forces [see Eqs. Q 
and (|2}]. This matter is not accounted for by the standard 
velocity-Verlet scheme. The DPD-VV, however, accounts for 
this complication in an approximate fashion by updating the 
dissipative forces [step (5) in TableU for a second time at the 
end of each integration step. This improves its performance 
considerably yet keeping it computationally efficient since the 
additional update of dissipative forces is not particularly time- 
consuming. In previous studies, the DPD-VV scheme has 
shown good overall performance 1 2lj| l22ll for which reason 
we have chosen it as the "minimal standard" to which other 
integrators are compared. 



B. Self-consistent velocity-Verlet integrator 

The update scheme of a self-consistent variant of DPD-VV 
is presented in Table ID This SC-VV algorithm llUllj de- 
termines the velocities and dissipative forces self-consistently 
through functional iteration, and the convergence of the iter- 



ation process is monitored by the instantaneous temperature 
kgT. This approach is similar in spirit to the self-consistent 
leap-frog scheme introduced recently by Pagonabarraga et al. 
1 20], which is the only other published self-consistent DPD 
integration scheme in addition to SC-VV (to the authors' 
knowledge). The SC-VV scheme has the advantage of be- 
ing very easy to implement as seen from Table U A recent 
study of the SC-VV scheme confirmed that it is a promising 
approach for interacting DPD systems |22]. That is partic- 
ularly the case for the structural properties using long time 
steps in dense systems. As a drawback, the SC-VV has no 
advantage in temperature control as compared to other meth- 
ods. For that, there exists a variant of the SC-VV integrator 
with a Nose-Hoover type additional thermostat i2lll22ll . 



C. Integrator by den Otter and Clarke 



(2) fi < — fi + ViAt 

(3) Calculate i^^{r,}, Fr{rj,Vj}, i^^{r,} 

(4) Calculate physical quantities 



TABLE II: The approach OC by den Otter and Clarke. Initialization: 
Calculate averages {Ff ■ Vi), {Ff ■ Ff) and (J^^ • J^^) either 
analytically or numerically. Then extract a and (3 from Eqs. ^6} and 
0, respectively. 



Parameter 


Model A 


Model B 


Model polymer system 




-7.528088 


-4.985245 


-7.566508 


{FF ■ fP) 


38.254933 


14.507340 


38.532832 


{Fi" ■ Fi") 


15.072610 


9.976453 


15.150489 



TABLE III: The values of parameters used in the OC integrator. 



In 2001 den Otter and Clarke |23] proposed an approach 
which uses a leap-frog algorithm with predefined variables a 
and (3. The idea in the OC-integrator, as it is called here, is 
to try to determine these parameters such that the effects due 
to the velocity dependence of dissipative forces are reduced 
as much as possible. The OC algorithm is described in 
Tableim in which the parameters a and (3 describe the relative 
weight of the random forces with respect to the dissipative 
and conservative ones. They are determined prior to the ac- 
tual DPD simulation by calculating the averages {Ff ■ Vi), 
{Ff ■ Ff), and {Ff^ ■ F[^) from an ensemble in which both 
the kinetic and configurational temperature equal the desired 
temperature ksT* |23]. Once they have been calculated, one 
obtains a and (3 with a desired time step At from the equa- 
tions 



1 



GAt 



(1 



-GAt 



with G = -(^;^-t/,)/(t;,-w,), and 



2ma{Fl^ ■ v,) - a^At {Fl^ ■ F^ 



(6) 



(7) 



Note that both a and (3 depend explicitly on At. Since the 



averages {F^^ 



D 



{Ff ■ Ff), and {Ff ■ Ff) can be derived 



analytically only for a limited number of systems, one usually 
has to calculate them from simulation (with a very small time 
step). This might be a problem in cases where properties such 
as density and temperature are varied over a wide range, since 
the parameters a and f3 should (at least in principle) be calcu- 
lated separately for all different conditions. Nevertheless, den 
Otter and Clarke have shown |23] that the OC algorithm per- 
forms well in both the ideal gas and a softly interacting DPD 
fluid. 

For the three models considered in this work, we deter- 
mined the parameters {Ff ■ v,), {Ff ■ Ff), and {Ff ■ Ff) 
separately for all cases. Their values are shown in Table IIIII 

Note that the parameters for the model polymer system in 
Tablelllllhave been determined by averaging over all particles 



in the system. An alternate way would be to find {M + I) 
different values for a and f3 by averaging over the solvent and 
monomer particles separately. However, this would require a 
major computational study prior to actual simulations and is 
therefore not feasible. Besides, it might be against the usual 
spirit as integration schemes are typically based on prefactors 
that are identical for all particles in a system. 



D. Shardlow's splitting method 

The most recent addition to DPD integrators has been in- 
troduced by Shardlow |24|. He applied ideas commonly used 
in solving differential equations to the case of integrating the 
equations of motion in DPD. The key idea is to factorize the 
integration process such that the conservative forces are cal- 
culated separately from the dissipative and random terms. Af- 
ter this splitting the conservative part can be solved using tra- 
ditional molecular dynamics methods, while the fluctuation- 
dissipation part is solved separately as a stochastic differential 
(Langevin) equation. To this end, Shardlow suggested two in- 
tegrators, called S 1 and S2, based on splitting the equations of 
motion up to first and second order, respectively. 

The formal approach involves a first order splitting using 
the Trotter expansion |4, 30] (integrator SI) and a second or- 
der splitting using the Strang expansion ]31] (integrator S2). 
The mathematical details of the derivation can be found in 
Sharlow's original article ]24]. It is important to notice that 
the power of the Trotter (Strang) expansion lies in the fact 
that it provides a general method for deriving symplectic algo- 
rithms. Importantly, the method works for both Hamiltonian 
and non-Hamiltonian systems; see Ref. L30.1 for a detailed dis- 
cussion of the Trotter expansion and its appUcations. 

Based on our simulation studies and the results presented 
in Ref. ]24] for a system related to the model B in the present 
work, the performance of both of the two splitting methods 
was found to be excellent with S2 displaying slightly better 
overall characteristics. Since the first order method (SI) is 
more efficient, we have chosen it to be on the spotlight. Here 
we use the same naming convention and call it S 1 . The algo- 
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(1) For all pairs of particles for which nj < rc 

(i) Vi< — Vi- |^7^^^(''ii)(«ij • eij)eijAt + ^:^auj{r,j)^ijeijVM 
1 1 



(ii) Vj 

(iii) Vi 

(iv) Vj 



2 



Vi + 



-auj{rij)^i.jeijVAi ' 



v, + lJ^F^-At 



(2) ^ - i 1 1 

(3) fi < — fi + ViAt 

(4) Calculate J;'^{f J } 

(5) lii^^c^^ 

(6) Calculate physical quantities 



1 j_ 

2 m 

2 m l+7ii)^(ri -)At 



JJ^>-J U V '2m l + 7tij^(rij)At 



At 
'At 



TABLE IV: The approach SI by Shardlow. 



rithm is presented in Table II Vl To the best of our knowledge, 
further tests of S 1 have not been reported yet. 

The curious fact that the algorithm (Table II V> appears 
asymmetric for particles i and j is a result of the use of the 
fluctuation-dissipation theorem and Newton's third law. The 
form in which the algorithm is presented keeps the notation 
otherwise symmetric. 



E. Lowe's approach - Lowe-Andersen method 

The approach introduced by Lowe in 1999 1 25] is presented 
in TablelVl which shows how one first integrates the Newton's 
equations of motion with a time step At, and then thermal- 
izes the system as follows. For all pairs of particles for which 
rij < rc, one decides with a probability FAt whether to take 
a new relative velocity from a Maxwell distribution. For each 
pair of particles whose velocities are to be thermalized, one 
works on the component of the velocity parallel to the line of 
centers and generates a relative velocity v^j ■ from a dis- 
tribution ^Ij''' ^2kBT* jm. Here are Gaussian random 
numbers with zero mean and unit variance. This approach 
has its origin in the Andersen thermostat |26], hence the name 
Lowe-Andersen method. 

The key factor in Lowe's method is the parameter l/F 
which describes the decay time for relative velocities. Since 
the condition < FAi < 1 is obvious, one finds that for 
Y At — 1 the particle velocities are thermalized at every time 
step, while for FAi w the model system is only weakly 
coupled to the thermostat. Thus the dynamical properties of 
the system can be tuned by the choice of F as shown by Lowe 
iQ. 

Although the present version of the algorithm follows the 
original one |25| and is based on the velocity- Verlet scheme, it 
is clear that other approaches such as the leap-frog are equally 
useful, if desired. Further, based on the work by Lowe |25], 
this approach seems very promising although it has received 
only little attention by far |32|. 

In the present work for the three model systems considered 
here, we set F such that the tracer diffusion properties of the 
fluid are similar with those of DPD systems with chosen cr in 
the limit At 0. In this fashion, we end up to a value of 



F = 0.745 for model A and to a value of F = 0.44 for model 
B. In the model polymer system we used the same value as in 
model A since the solvent is described in a similar fashion in 
both cases. 



(1) v^'^v, + \^F? At 

(2) Ti < — fi + ViAt 

(3) Calculate f;'='{r J } 

(4) v^^v^ + ll^FfAt 

(5) For all pairs of particles for which r,:j < 

(i) Generate v°j ■ e,;j from a distribution ^'j' \j2kBT* jm 

(ii) 2Aij = eij{vij - Vij) ■ Sij 

(iii) Vi < — Vi + Aij 

(iv) Vj < — Vj — Aij 
with probability FAt 

(6) Calculate physical quantities 

TABLE V: The approach by Lowe using Gaussian distributed ran- 
dom numbers ^^^^ from a distribution ^'j' ^2kBT* jra. 



IV. PERFORMANCE OF INTEGRATORS 

A. Physical quantities studied 

We characterize the integrators by studying a number of 
physical observables. After equilibrating the system, we first 
calculate the average kinetic temperature 

whose conservation is one of the main conditions for reliable 
simulations in the canonical ensemble. 

Next, in the cases of models A and B, we examine the ra- 
dial distribution function g[r) |33] which is one of the most 
central observables in studies of structural properties of liq- 
uids and solids. For the ideal gas (model A), the radial dis- 
tribution function provides an excellent test for the integrators 
since then g(r) = 1 at the continuum limit. Therefore, any 
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deviation from unity has to be interpreted as an artifact due to 
the integration scheme employed. 

In model B, in which conservative interactions are present, 
there are no exact theoretical predictions for g(r) that would 
allow a straightforward comparison of different integration 
schemes. A comparison is possible, though, in terms of phys- 
ical observables such as the compressibility and the coordina- 
tion number that are based on integrating g{r). In the present 
study, we have chosen to consider the coordination number 
defined as 



drg(r) 



(9) 



where p is the particle number density of the system and ri 
is the radial distance at which g{r) has its first minimum after 
the leading (first) peak. 

The radial distribution function reflects equilibrium (time- 
independent) properties of the system. To complement the 
comparison of different integrators, we also consider the 
tracer diffusion coefficient 



lim ^Am) 



(10) 



which can provide us with information of possible problems 
on the dynamics of the system. Here fi (t) is the position of a 
tagged particle in models A and B, and the mean-square dis- 
placement is then averaged over all particles in a system to get 
better statistics for Dt- In the model polymer system, ri{t) 
describes the center-of-mass position of the polymer chain via 



1 



(11) 



where the index runs over monomers in a chain. 

For the model polymer chain, we further calculate the ra- 
dius of gyration Rg = w (R^) defined as 



1 



M 



M ^ 

i=l 



(12) 



which shows that Rg is a measure of polymer size. It is actu- 
ally one of the most central quantities in polymer science and 
therefore serves as an excellent measure for our purposes. 

In the following, the errors are stated in the figure captions 
and given as the magnitude of standard deviation. 



B. Results for model A 

As discussed above, model A is characterized by the ab- 
sence of conservative forces, and thus any artifacts arising 
from the velocity-dependent forces are expected to be pro- 
nounced in this model. To study this possibility, we first dis- 
cuss the deviations of the observed kinetic temperature (fc^T) 
from the desired temperature ksT* . The results for {ksT) 



1.00 



-e-e — e-e- 



-k — i -I 




0.01 



0.1 



^t 



FIG. 1: Results for the deviations of the observed temperature 
(kBT) from the desired temperature ksT* = 1 vs. the size of the 
time step At in model A. The error in (fcsT) is of the order of 10^''. 



shown in Fig. [J indicate that DPD-VV and SC-VV are rea- 
sonably good at small time increments, but larger time steps 
lead to major deviations from the desired temperature. For 
OC, (ksT) decreases monotonically with At for At < 0.15, 
after which the temperature increases rapidly. Nevertheless, 
the deviation is greater than in the case of DPD-VV but con- 
siderably smaller than for SC-VV. The Shardlow S 1 integra- 
tor, however, has very good temperature control and the de- 
viations remain less than 0.5 % up to At — 0.2. The best 
temperature control is found for the method by Lowe, how- 
ever, yielding (ksT) = ksT* for all time steps At. In this 
case, we have extended the studies further and tested the be- 
havior of (ksT) with various values of F between 0. 1 and 10, 
but the conclusions remain the same. 

Results for g(r) are shown in Fig.|2] We find that the de- 
viations from the ideal gas limit g{r) ~ 1 are pronounced 
for OC, indicating that even for small time steps this integra- 
tion scheme gives rise to unphysical correlations. The perfor- 
mance of DPD-VV is considerably better, although artificial 
structures are yet rather pronounced, while SC-VV and SI 
lead to a radial distribution function that is close to the theo- 
retically predicted one. Completely structureless g{r) is found 
only for the integrator by Lowe, however. Again, in this case, 
we have tested the behavior of g{r) with various values of F, 
but the results remain the same. This confirms the expecta- 
tion that F does not influence the equilibrium properties of the 
system. 

The results for the diffusion coefficient Dt in Fig.|3]are es- 
sentially consistent with the conclusions above. The integrator 
OC is not very useful in a system of the present kind, since it 
seems to lead to substantial deviations from the expected be- 
havior. The SC-VV and the integrator by Lowe perform much 
better, while the integrators S 1 and DPD-VV are most stable 
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FIG. 2: Radial distribution functions g{r) with several values of time step At in model A: (a) At = 0.01, (b) At — 0.05, and (c) At = 0.1. 
The error in g{r)is greatest at r = 0.01, where it takes the value of 0.01. 




0.01 0.1 0.01 0.1 

At At 



FIG. 3: Results for the tracer diffusion coefficient 

Dt { At) /DriO .01) vs. the time step At in model A. The 
error in Dt (At) / Dt (0.01) is of the order of 0.001. 



FIG. 4: Results for the deviations of the observed temperature 
{ksT) from the desired temperature ksT* = 1 vs. the size of the 
time step At in model B. The error in (ksT) is of the order of lO"**. 



in this case. 



C. Results for model B 

Results shown in Fig. |4] for the observed kinetic tempera- 
ture {ksT) reveal that the differences between the integration 
schemes deviations are weaker in model B than in model A. 
As it turns out below, this conclusion is generic and holds for 
all quantities studied here. 

The deviations of {ksT) from the desired temperature 
ksT* are very minor for all integrators at small time steps 



At < 0.01. Differences between the integrators become ev- 
ident only at larger time steps. We first find how the temper- 
ature in SC-VV first decreases, then has a pronounced min- 
imum around At w 0.05, after which (fcsT) increases very 
rapidly and the integrator eventually becomes unstable. The 
temperature conservation of OC is also relatively poor at time 
steps beyond At — 0.025. Best performance in this respect 
is found for the remaining integration schemes DPD-VV, SI, 
and Lowe, whose behavior is quite similar and comparable to 
each other. 

Demonstrative results for the radial distribution functions 
in model B are shown in Fig.|5] It is clear that large time steps 
lead to major problems with regard to pair correlations. This 
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FIG. 5: Radial distribution functions g(r) witli time steps At — 0.01 (on tlie left) and At = 0.1 (on the right) in model B. In addition to the 
full curves, two sets of the same data on an expanded scale are also given to clarify the deviations between different integration schemes. The 
error in g{r)is greatest at r = 0.01, where it takes the value of 0.001. 
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is particularly clear at small distances (r < 0.2). To quantify 
these changes, we calculated the coordination number [see 
Eq ©]- The results shown in Table [Vl] highlight the fact that 
all integration schemes converge to the same result at small 
time steps At < 0.05, while for larger time steps there are sig- 
nificant deviations from the correct behavior found in the limit 
At — > 0. However, it is somewhat surprising that the coordi- 
nation numbers obtained by different integration schemes are 
essentially similar within error limits, while based on g(r)'s 
there are noticeable differences between the pair correlation 
properties of different integrators. This is a consequence of 
the fact that due to the integration in Eq. (|9}, the deviations in 
g(r) toward too small and too large values compensate each 
other A similar effect has been observed recently in com- 
pressibility 1 22], which is also defined as an integral over 5 (r). 

Tracer diffusion data shown in Fig.|6lis consistent with the 
conclusions above. For small time steps the results of all in- 
tegration schemes are comparable, while for large time steps 
we can find how the differences become more and more pro- 
nounced. The scatter in the data does not allow us to make 
conclusive statements of the relative merits of the integrators, 
however. 




0.01 



0.1 



At 



FIG. 6: Results for the tracer diffusion coefficient 

Dt (At) /Dt (0.005) vs. the time step At in model B. The 
error in Dt (At) / Dt (0.005) is of the order of 0.001. 



D. Results for the model polymer system 

Before we consider the results for the model polymer sys- 
tem, we would like to emphasize certain similarities it has 
with model A. Namely, in both cases the solvent is an ideal 
gas governed by dissipative and random forces only. Further, 
the model polymer system is dilute (more than 99.5% of the 
particles in a system are solvent particles) and there are no 
conservative interactions between the monomers and the sol- 
vent particles. This suggests that the artifacts due to integra- 
tion schemes in the model polymer system would be essen- 
tially similar to those found in model A. It turns out below, 
however, that this is not the case. 

Results shown in Fig. 0a) indicate that the observed ki- 
netic temperature (fcsT) only rather weakly depends on the 
integration scheme. Results for SI, Lowe, and DPD-VV are 
all within one percent up to At « 0.02 above which the inte- 
gration schemes become unstable. The OC integrator is some- 
what less reliable in this case, as it leads to a monotonous in- 
crease of (ksT), thus differing rather clearly from the results 
of other integration schemes. 

A comparison of Figs.^and^Ia) provides one with an in- 
triguing view of effects that arise from the hybrid approach. 
We first note that Lowe's approach as well as SI are equally 
good, in agreement with the conclusions made in model A. 
However, in model A the integration schemes were found to 
be stable up to very large time steps on the order of At k, 0.4, 
while in the model polymer system the largest time steps pos- 
sible are about 0.02. This is due to the hard conservative 
monomer-monomer interactions used in describing the poly- 
mer chain, for which reason the size of the time step has to be 
reduced considerably as compared to model A. This suggests 
that, for practical purposes in hybrid models, one should se- 
riously consider integration schemes with two different time 



steps, one for the solvent and another for the polymer degrees 
of freedom. Another interesting feature concerns the behavior 
of (fc^T) in the case of OC. In model A, the observed ki- 
netic temperature decreased monotonously down to timesteps 
At « 0.15, while in the model polymer system the trend is the 
opposite. This is consistent with our results for model B, and 
suggests that the artifacts due to integration schemes depend 
significantly on the interactions chosen for the system. 

To gain further insight into the performance of the integra- 
tion schemes in a hybrid approach, we studied a number of 
physical quantities that specifically characterize the properties 
of the polymer chain. First, we studied the observed kinetic 
temperature of the polymer, defined as 

(A:i3r)chain= ^^^t^^y (13) 

This quantity characterizes the thermal fluctuations of the 
polymer chain, and therefore if there are any serious prob- 
lems due to the integration schemes, then we expect that they 
are manifested in the behavior of ( fc^r) chain- 

We find that the results in Fig. 0b) are wholly consistent 
with those presented in Fig. 0a). Essentially, this implies that 
the temperature deviations in the whole system arise from the 
hard interparticle interactions used to describe the solute even 
though the system is dilute. We think that this finding is of 
generic nature and applies to both hybrid models and other 
DPD simulations in which all interactions are approximately 
of equal magnitude. In particular, it allows us to suggest that 
the artifacts due to integration schemes in DPD are predomi- 
nated by the interactions that dictate the size of the time step. 

In systems with hard conservative interactions the time step 
must be small. Otherwise, gradients in forces become too 
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Integrator At = 0.005 At = 0.010 At = 0.050 At = 0.075 At = 0.100 



DPD-VV 


25.36 


25.37 


25.74 


26.72 


28.98 


SC-VV 


25.43 


25.54 


25.62 


26.78 


28.57 


oc 


25.46 


25.51 


25.71 


26.69 


28.74 


SI 


25.53 


25.41 


25.73 


26.71 


28.44 


Lowe 


25.35 


25.42 


25.59 


26.95 


29.12 



TABLE VL Results for the coordination number Nc in model B. Error bars are about ±0.05. 
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FIG. 7: Results for the deviations of the observed temperature (kBT) from the desired temperature ksT* = 1 vs. the size of the time step At 
in the model polymer system. In (a) we show the observed temperature of the whole system, while the results in (b) correspond to the polymer 
chain only (see text for details). For the whole system, the error is of the order of 10"'*, while for the polymer chain the error is 0.004. 



large and the system becomes unstable. Consequently, if the 

conservative force is stronger than the dissipative and random 
forces the artifacts due to integration schemes are driven by 
the conservative forces just like in classical molecular dynam- 
ics simulations. Naturally, the velocity-dependent dissipative 
forces are still playing a role but their effect is not as important 
as the influence of conservative interactions. From a practical 
point of view, this means that there is no particular reason to 
use an integrator which accounts for the velocity dependence 
of dissipative forces. 

On the other hand, if all interactions are soft, or if the con- 
servative forces are weak compared to the dissipative forces, 
then it is plausible that the velocity-dependence of dissipative 
forces is the underlying reason for artifacts due to the inte- 
gration procedure. This is the case when the time step is de- 
termined by the dissipative force rather than the hard-core of 
the conservative potential. In this situation the quality of the 
integration scheme is very important, and the velocity depen- 
dence of the dissipative forces has to be accounted for by the 
integration scheme. It is clear that this matter warrants atten- 



tion and should be accounted for in all subsequent studies by 
DPD. 

Further studies of the radius of gyration and the diffusion 
coefficient of the polymer chain revealed the expected and 
undesired fact that computational studies of a dilute polymer 
system in an explicit solvent are very time consuming, and 
therefore the error bars remained rather large despite major 
computational efforts. Consequently, we found that the re- 
sults for Rg and Dt of different integrators were essentially 
equal within error limits (data not shown). It is likely that 
more extensive calculations would have expressed deviations 
between different integration schemes, but we concluded that 
such studies were not worthwhile. 



E. Computational efficiency 

Besides the strength of the artifacts, we have paid attention 

to the computational efficiency of the integration schemes by 
calculating the cpu time needed for a single time step At. Al- 
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Integrator Cpu time (seconds) 

DPD-VV 0.0363 ± 0.0005 

SC-VV 0.1010 ±0.0009 

OC 0.0251 ± 0.0005 

SI 0.0256 ± 0.0005 

Lowe 0.0143 ± 0.0005 

Verlet list 0.0293 ± 0.0003 



TABLE VII: Results for the computational efficiency of the integra- 
tion schemes. Shown here are results for integrating the equations of 
motion over one time step of size At = 0.05, although the results 
have been averaged over 1 000 consecutive steps. For the purpose of 
comparison, the time needed to update the Verlet neighbor list has 
also been given; the time shown here corresponds to its minimum 
value when the list is small since it is updated after every time step. 
(Simulation parameters: fcsT = 1, A — 25, a — 3, and p = 4.) 



though this depends on various practical matters such as the 
computer architecture, the implementation of the algorithms, 
and the size of Verlet neighbor tables, we think this approach 
can provide one with the essential information of the relative 
speed of the integrators tested in this work. 

The tests for efficiency have been carried out on a Compaq 
Alpha Station XPIOOO with a 667 MHz processor. To this end, 
we used model B (with 4 000 particles at a density of p = 4) 
with standard Verlet neighbor tables ll28il and a time step of 
At = 0.05. 

We calculated the cpu time needed for the integration of a 
single time step (averaged over 1 000 consecutive steps). The 
times needed to update the Verlet neighbor tables (or to cal- 
culate any physical quantities) are not included in these re- 
sults. Thus the DPD-VV and SC-VV schemes were consid- 
ered over steps (l)-(5) in Table|I] the OC approach over steps 
(l)-(3) in Table HII Shardlow's approach over steps (l)-(5) 
in Table Hy] and Lowe's integration scheme over steps (1)- 
(5) in Table|Vl In the case of SC-VV, steps (4b) and (5) were 
repeated six times which guaranteed self-consistency in this 
case. Note that Lowe's approach is based on normally dis- 
tributed random numbers while other integration schemes use 
uniformly distributed ones. Since there are various methods 
available for the generation of normally distributed random 
numbers, this may affect the efficiency of Lowe's approach to 
some extent |34]. 

The results shown in Table lVlll indicate that Lowe's method 
is substantially faster than OC and SI, which in turn are 
clearly faster than DPD-VV. Finally, the result that SC-VV 
is considerably slower than DPD-VV is not surprising due to 
the iteration process for the velocities and dissipative forces. 

When these times are compared to each other, one should 
also keep in mind that DPD-VV, SI, and Lowe's method are 
significantly easier to deal with compared to SC-VV and OC. 
In SC-VV, one needs to find the sufficient number of itera- 
tions prior to actual simulations to guarantee self-consistency. 
In OC, the preliminary work required prior to simulations is 
even more extensive, as one has to determine the parameters 
a and (3 for a given system under desired thermodynamic con- 
ditions. This task may indeed take some time. 



V. DISCUSSION AND CONCLUSIONS 

In this work, we have tested several novel schemes on 
an equal footing through DPD simulations of three different 
model systems. The first of the models coiTesponds to a case 
where conservative interactions play no role, while the sec- 
ond model describes fluid-like systems with relatively strong 
but soft conservative potentials. Finally, the third model aims 
to characterize the quality of integration schemes in a hybrid 
approach for a dilute polymer system. 

Of the integration schemes considered here, DPD-VV and 
SC-VV have recently been examined in Refs. i2lll23l . The 
results of the present study are consistent with previous find- 
ings: DPD-VV exhibits good overall performance, indicat- 
ing that it presents a relatively accurate means to integrate the 
equations of motion at a reasonable computational cost. 

Of the previously untested methods the Otter-Clarke (OC) 
method |23| is fast, performing especially well in interacting 
systems in which conservative forces are important. A draw- 
back is that the parameters a and (3 need to be determined 
prior to actual simulations through time-consuming precur- 
sory simulations with a very small time step. We note, how- 
ever, that the properties of the OC scheme are in fact relatively 
insensitive to slight changes in a and f3. Thus, for example, 
studies of the model polymer system using specifically deter- 
mined a and f3 values yielded results almost identical with 
studies based on the parameters of model A. 

The Shardlow SI integrator |24] is possibly the brightest 
star in this work. It performed very well in all models, and it 
is fast and rather easy to implement. We feel that it presents 
the best choice of integration schemes within the "usual" con- 
ceptual framework of DPD. 

Interestingly, however, we have also found that the elegant 
and conceptually distinct method of Lowe |25] performed ex- 
cellently and is easy to implement. Furthermore, and what 
is important when Lowe's method is compared to Shardlow's 
integration scheme, it provides an alternative and a very at- 
tractive description of dissipative particle dynamics. Thus, a 
direct comparison of SI and Lowe's method is not meaning- 
ful. Instead, we discuss the pros and cons of these two ap- 
proaches. 

The usual DPD description is based on the idea that soft 
matter systems can be described in terms of softly interacting 
particles with some of the degrees of freedom coarse grained 
out and replaced with random noise coupled to dissipation. 
Temperature conservation is achieved through the fluctuation- 
dissipation theorem and the correct hydrodynamic behavior is 
guaranteed by momentum conservation 1 10]. Various studies 
have extended these ideas further. For example, Flekk0y et 
al. developed a DPD framework starting from a microscopic 
description Issl l36ll . Espaiiol and coworkers, in turn, stud- 
ied the dependence of transport properties of DPD fluids on 
the length and time scales |37] and a generalization of DPD 
to energy conserving systems |38]. DPD has recently been 
used together with molecular dynamics to coarse grain aque- 
ous salt solutions |39] in which the effective interactions used 
in DPD simulations were obtained from MD simulations by 
the inverse Monte Carlo procedure |40]. The last ten years 
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FIG. 8: a) Velocity autocorrelation function for Model A in Lowe's method. The error is of the order of 10 b) Velocity autocorrelation 
function for Model B in Lowe's method. The error is of the order of 10"*. 
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FIG. 9: Velocity autocorrelation function for Model B in Lowe's method with the product TAt fixed to (a) 0.005 and (b) 0.025. The error is 
of the order of 10"". 



have been very successful on both the analytical and the com- 
putational fronts-the theoretical basis of DPD is now well es- 
tablished, and the number of applications has increased at a 
steady pace. 

Lowe's |25] approach is a very recent inception and thus far 
has received limited attention. Although the theoretical foun- 
dations of Lowe's method have yet to be fully worked out, it 
offers promising aspects that are not obvious in the traditional 
DPD description. To clarify these aspects, let us first remind 
ourselves that Lowe's method does not include dissipation in 
the usual sense. Rather, it is based on a thermostat that ther- 
malizes the velocities of pairs of particles at a rate which de- 
pends on the dynamical parameter F. This parameter tunes 
the dynamical properties of the system. Lowe pointed out that 
the soft interactions used in DPD lead to a situation where the 
ratio of the kinematic viscosity and the diffusion coefficient 
of solvent particles (known as the Schmidt number Sc) is of 
the order of one. This value corresponds to a situation often 
found in gases, while in fluids Sc ~ 10^ or even larger. To get 



closer to more realistic values for Sc, one can reduce the dif- 
fusion rate by using harder interparticle interactions, but this 
is against the philosophy of DPD and would reduce some of 
the benefits of the DPD approach. 

Lowe's approach is very different in this respect. It allows 
one to adjust the viscosity of the system to a desired value by 
varying the dynamical parameter F while the diffusive prop- 
erties are not considerably affected since the conservative in- 
terparticle interactions remain soft. As a result, the Schmidt 
number can obtain values as large as 10' When com- 

pared to the usual DPD description, this implies that Lowe's 
approach may be more feasible for describing hydrodynamic 
systems in which one needs to worry about the time scales 
of momentum diffusion and mass transfer with respect to the 
size of the colloidal particle. 

There still remains the issue of the practical viability of 
Lowe's approach, since we are not aware of any applica- 
tions where the method by Lowe has been used. However, 
we are positive that this approach is a promising technique. 
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For example, we have recently applied Lowe's method to mi- 
crophase separation of block copolymers in the spirit of Groot 
and Madden |41], and it turned out that Lowe's method was 
able to reproduce their results. Finally, in Figs. |8]and|9]we 
show how the velocity autocorrelation function depends on 
the choice of F for models A and B (Fig.|8|i, and how it is af- 
fected when r is varied but keeping the product FAt constant 
in the case of model B (Fig.|9jl. As discussed above, it is clear 
that large values of F lead to faster decay. However, the qual- 
itative behavior of the velocity autocorrelation function is not 
seemingly affected, as illustrated by Fig.|9] and the effect of F 
on the diffusion coefficient was found not to be important for 
the studied combinations of F and At. 

To conclude, we have studied the performance of various 
novel integration schemes that have been designed specifi- 
cally for DPD simulations. We have tested these integration 
schemes in three different model systems by varying the na- 
ture of interactions and found that the artifacts due to the inte- 
gration scheme are essentially driven by the interactions that 
dictate the size of the time step. Thus, the artifacts and the 
performance of integrators are model dependent. Overall, we 
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